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In this paper, we present a computationally efficient method for including fluid-solid interactions 
into direct numerical simulations of the Navier-Stokes equations. This method is found to be as 
powerful as our earlier formulation [J. Comp. Phys., vol. 249: 243 (2015)], while outperforming the 
earlier method in terms of computational efficiency. The performance and efficacy of the presented 
method are demonstrated by computing contact angles of droplets at equilibrium. Furthermore, we 
study the instability of films due to destabilizing fluid-solid interactions, and discuss the influence 
of contact angle and inertial effects on film breakup. In particular, direct simulation results show 
an increase in the final characteristic length scales when compared to the predictions of a linear 
stability analysis, suggesting significant influence of nonlinear effects. Our results also show that 
emerging length scales differ, depending on a number of physical dimensions considered. 



I. INTRODUCTION 


Stability of thin films, in particular at the nanoscale, is relevant in a variety of applications where film 
breakup and dewetting are important. One example of such applications involves harnessing instabilities 
of nanoscale liquid metal films in order to create arrays of nanoparticles via self- or directed-assembly [1- 
4]. These methods involve the deposition of thin, solid metal films, that are subsequently exposed to laser 
pulses; breakup leads to the formation of drops, which then solidify to form nanoparticles. The applications 
of such nanoparticles are wide, ranging from solar cells to liquid crystal displays, among others [5-10]. 
Such nanofilm breakup presents numerous computational challenges that may result from complex initial 
geometries [11, 12]; however, rupture of a simple flat film due to fluid-solid interactions also demonstrates 
rich behavior that has been studied extensively in the past (see, for example, [13-15], or [16] for a review). 

To model the film breakup and the consequent dewetting phenomena, it is necessary to include a desta¬ 
bilizing mechanism: if such a mechanism is not included, a continuous film does not break up. In this 
context, the long-wave formulation (L-W) is usually considered, since it simplifies the underlying mathemat¬ 
ical model, reduces dimensionality of the problem, and comes at significantly reduced computational cost. 
Liquid-solid interaction forces are often included in the L-W as conjoining-disjoining pressure, leading to a 
prewetted (often called ‘precursor’) layer in nominally ‘dry’ regions. This approach effectively removes the 
‘true’ contact line, consequently avoiding the non-integrable shear-stress singularity at the moving contact 
line (see e.g. [17-19]). Another approach to alleviate the moving contact line singularity is to relax the no-slip 
condition and instead assume the presence of slip at the liquid-solid interface (see e.g. [20-23]). Slip at the 
solid surface for fluid-fluid systems has also been studied in the context of molecular dynamics simulations 
(see e.g. [24-27]). Both slip and disjoining pressure approaches have been extensively used to model a variety 
of problems including wetting, dewetting, and film breakup, in particular in the context of polymer [28, 29] 
and metal [1, 3, 4, 30-33] Aims. 

The L-W, however, has its limitations, in particular regarding the requirements of small interfacial slopes 
and negligible inertial effects. For example, the L-W may not be sufficient to provide quantitatively precise 
predictions regarding instability development for metal Aims at nanoscales, where contact angles are large, 
and inertial effects considerable [34-37]. Regarding the small slope requirement, it is often argued that the 
results of the L-W are reasonably accurate even if this requirement is not strictly satisfied; in addition, 
various improvements of the L-W are available, see, e.g. [38] . Still, since it is not always possible to compare 
results of simulations directly to experiments, it is difficult to judge to which degree the L-W based results 
quantitatively describe the process of thin Aim breakup. Regarding inertial effects, although the L-W can 
be extended to include them (see [39, 40] for reviews of this topic), the resulting formulations are not 
straightforward, and are limited in the range of Reynolds numbers that can be considered. 

To be able to accurately model the problems where the assumption of small slopes is not satisfied, or 
where inertial effects are important, it is necessary to go beyond the L-W, and consider direct solutions of 
Navier-Stokes equations, that also include fluid-solid interaction forces. It is therefore essential to design 
a computational method to include the following: (i) precise resolution of long and short range fluid-solid 
interactions; (ii) allow non-negligible inertial effects and large contact angles; (hi) provide spatially and 
temporally converged solutions with a reasonable computational cost. 

In our earlier work [41], we developed a model for inclusion of fluid-solid interaction forces in a Navier- 
Stokes solver. However, that model requires significant computational effort, and it is not practical for use 
in more complex scenarios, such as three dimensional (3D) Aim breakup. In the current paper, we present 
a computational approach that is significantly more efficient. The present approach permits the analysis of 
complex evolution of the rupturing and dewetting process. Our results reveal that the Aim initially evolves in 
accordance to the predictions of the linear stability analysis (LSA) based on the L-W, but departs from the 
LSA prediction beyond the onset of the instability. We also discuss the influence of the number of physical 
dimensions considered, by showing differences in the rupture process between two and three dimensions 
(2D and 3D). Despite the fact that a variety of other computational methods have been considered in the 
context of wetting/dewetting (see e.g. [42-44]), to our knowledge, the numerical scheme presented in this 
work provides the only available framework that satisfies all three requirements listed above to quantitatively 
describe the dynamical evolution of thin Aim instability including rupture. 

This paper is organized as follows. We begin in Sec. II by giving an overview of the computational 
framework and a brief review of the L-W. We also briefly discuss the LSA of the L-W equations of an 
initially flat Aim on a substrate. In Sec. Ill, we compare the results for the computed contact angle with 
available results. Then, we present the results of the Aim instability in linear as well as nonlinear regimes 
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in 2D. We study the underlying physical process, concentrating in particular on the role of the inertia and 
contact angle to better understand the regimes of the validity of the LSA. Finally, we present simulations of 
nonlinear evolution and breakup of 3D films and contrast the results with 2D ones. The main conclusions 
are summarized in Sec. IV. 


II. MODEL 
A. Derivation 

We consider a solid phase occupying a half-infinite region ^ < 0, above which there is a region occupied 
by two fluids that we conventionally refer to as a vapor phase (variables subscripted v) and a liquid phase 
(variables subscripted /). The words ‘vapor’ and ‘liquid’ are used purely conventionally, and do not necessarily 
signify anything about the physical properties, except that the two phases may differ in terms of their densities 
and viscosities. 

Each particle of fluid phase i interacts with the solid substrate by means of a Lennard-Jones type poten¬ 
tial [45] 


where r is the distance between the two particles, and is the scale of the potential, such that Eq. (1) 
has a minimum at r = We can obtain the total potential energy in phase i due to this 

interaction by the following expression 

/ OO ^0 nOO 

/ / (pisrisdxdydz, ( 2 ) 

OO j —OO j —OO 

where rii and are the densities of particles in fluid phase i and the solid substrate, respectively. Performing 
the integration, we obtain 


where 
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Equation (3) gives the total potential per unit volume in fluid phase i due to the interaction with the solid 
substrate. The quantity h* is conventionally referred to as the equilibrium film thickness in the literature [18, 
46], and we will continue using this convention here. The term equilibrium film thickness arises because h* 
is the thickness of a film that represents an energetic minimum due to Eq. (3), and in this model completely 
wets the surface of the substrate. 

We consider the inclusion of the potential in Eq. (3) in the Navier-Stokes equations for two phases. Eor 
clarity, we will refer to two phases as the liquid phase (subscript /), and the vapor phase (subscript v)^ 
i.e. i = although the present formulation applies to any two fluids. In order to do this, we introduce a 
characteristic function x(x, z)^ which takes the value of 1 inside of the liquid phase, and 0 inside the vapor 
phase. The interface between these two phases is assumed to be sharp, so that y changes discontinuously at 
the interface. The governing equations consequently become 

Pix)^ = -Vp + V • [p(x) (Vu + Vu^)] + jKSsn - /C(x)VF(y), ( 6 ) 




m=p — 3, n = q — 3. 


( 5 ) 
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V-u = 0. 


( 7 ) 


Here, p is the (phase dependent) density, p is the viscosity (phase dependent), p is the pressure, and u is 
the velocity vector. Also, D/Dt = dt is the material derivative. The surface tension is included as 

a singular body force [47]. Here 7 is the coefficient of surface tension, hz is the interfacial curvature, 6s is a 
delta function centered on the interface, and n is a normal vector for the interface pointing out of the liquid. 
The quantities depend on y via the following 

P = PiX + Pv{i-x), 

P = PiX + Pvi^ -x), 

^ = ^IsX + ^«s(l “ x)- 


In order to nondimensionalize, Eqs. ( 6 ) - (7), we introduce the length scale L, and define the time scale 
as the capillary time, r = /i/T/y, since, for the parameter sets considered in this paper the dynamics is 
dominated by viscous and surface tension effects. The dimensionless variables are defined as follows 


X=-; y=-- z=-; 



h* 


h* ^ ur ^ Lp ^ 

—; u=—; p =—; n = Ln] p 

L L ^ 
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With these scales, and dropping the tildes, the dimensionless Navier-Stokes equations become 


(Oh^)-= -Vp + V • [p(Vu + Vu^)] + K(5,n - KyF{y), ( 8 ) 

where the Ohnesorge number is defined via Oh^ = p/{pi^L), with L chosen according to the problem under 
consideration, and 


K = Kigx + K,g{l -x) = ^ = ^ {ICux + - X)) ■ 

Equation (7) is invariant with this scaling. In [41], we solved Eq. ( 8 ) directly. In that paper, we described 
two methods that we here refer to as ‘body force’ methods; in these methods, we discretized the fluid- 
solid interaction term differently, but we found both methods to be substantially similar in their overall 
properties. Eor consistency, when we refer to the ‘Body Eorce’ (B-E) method in this paper, we exclusively 
refer to Method H in [41]. As described in our previous work, the fluid-solid term presents a challenge for 
body force methods, in that the potential is divergent as p —> 0 . Consequently, parts of the domain near the 
substrate require a high mesh resolution in order to obtain reasonably accurate results. This requirement 
significantly limits the use of adaptive meshes, that are crucial for the purpose of improving the performance 
of direct simulations. In 2D, simple problems are feasible, however, in 3D, the computational cost of resolving 
so much of the domain makes the simulations impractical. 

In this section, we show that the computational task is dramatically simplified by reformulating the body 
force term in Eq. ( 8 ) as a force which acts only on the interface. Eirst, we define 


p* =p + KigxF{y) + Kyg{l - x)F{y), 


SO that 


-Vp* = -Vp - {Kux + K,g{l - x))VF - {K,g - Kig)5gnF{y), 

where Sgii = — Vy, in a distributional sense. Substituting into Eq. ( 8 ), we obtain what we refer to as the 
‘Reduced Pressure’ (R-P) formulation 

(Oh^)-= -Vp* + V • [p(Vu + Vu^)] + (« + >cF{y))5gn, (9) 
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where k = {Kvs — Kis) (note that only this difference is relevant, rather than individual values of Kys and 
Kis). 

Simple energy arguments [41] can be used to show that the ffuid-solid interaction term in Eq. (9) gives 
rise to an equilibrium contact angle ^eq for the liquid phase given by 

^ _ (1 - cos6>eq) / (m - l)(n - 1) 
h* \ m — n 

The exact meaning of ^eq requires some clarification. For the ffuid-solid interaction considered (with conjoin¬ 
ing and disjoining components), there is a film of thickness h* that wets the entire substrate, with a smooth 
transition from the droplet to the equilibrium film. We measure contact angles by fitting a circular profile 
to the droplet profile ‘away’ from the transition region; the angle at which this circular profile intersects 
the equilibrium film is taken to be the contact angle. The angle ^eq is the equilibrium contact angle in the 
following sense: for a drop at equilibrium with a vanishingly small h*, applying the above fitting procedure 
yields 6>eq. For non-vanishing but small h*, a drop at equilibrium will have a (slightly) different contact 
angle, which we refer to as ^num- 

In addition to the contact angle, the ffuid-solid interaction gives rise to another phenomenon that we 
consider in this paper: the spontaneous rupture of a thin liquid film. Our account of film rupture comes 
from [18], where this phenomenon is studied in the context of the L-W. For completeness, in Sec. IIC, 
we outline the L-W based approach for inclusion of the ffuid-solid interaction using the disjoining pressure 
model. 



B. Computational Implementation 


The R-P method discussed in this paper possesses three important strengths. First, in contrast to con¬ 
ventional methods for implementing the contact angle, that impose it essentially as a constraint that the 
solution needs to satisfy, it possesses the most important feature of the B-F method discussed in [41]: it 
includes long range fluid-solid interaction, and therefore spontaneous film breakup can occur. This feature 
is crucial if one wants to analyze instability of fluid films on nanoscale. The second major advantage of the 
formulation presented in this work is that by absorbing the entire contribution of the ffuid-solid interaction 
into a surface force, the main weakness of the B-F method is avoided. As discussed previously [41], the 
B-F method significantly limits the usefulness of adaptive mesh refinement. Since the force term in the B-F 
method applies everywhere in the fluid domain, and is singular as ^ ^ 0, the entire domain near y = 0 
requires a finer mesh in order to resolve the force; a typical adaptive mesh for simulations from [41] using 
the B-F method is shown in Fig. 1(a). In the R-P method, the force due to the fluid-solid interaction is 
applied as a surface force; consequently, high resolution is only required at the interface. Figure 1(b) shows 
an adaptive mesh used for simulations using the R-P method in this paper. A further advantage of the 
R-P method is that it is relatively simple to implement, owing to the fact that it can be formulated as a 
modification of the curvature in the surface tension term in the Navier-Stokes equations. 

We solve Eqs. (9) using the software package Gerris [48], described in detail in [49]. The mesh consists of 
a quad tree (in 2D) or an octree (in 3D), that decomposes the domain into square control volumes, which 
we refer to as cells (see Fig. 1). We use an adaptive mesh for all simulations, refining the interface to a 
resolution of Amax, and to lower resolutions far away from the interface. 

The interface is tracked using the Volume of Fluid (VoF) method. The VoF method replaces the charac¬ 
teristic function y of Eq. (6) with its average over each cell. This cell average is referred to as the volume 
fraction, T, and denotes the fraction of each computational cell occupied by the liquid phase. The VoF 
method reconstructs the interface as a sharp interface which is piecewise linear in each cell. Our fluid-solid 
interaction is included as a modification of the curvature in the surface tension term described in [49]. In 
this method, the curvature, /^, is calculated at cell centers, and our method replaces it with k - 1 - KF{y), 
where y is the ^-coordinate of the midpoint of the interface in an interfacial cell. The boundary conditions 
are as follows: on the solid substrate, we impose no-slip and T = 1; on all other boundaries, we impose a 
homogeneous Neumann boundary condition on all variables. 
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FIG. 1: (a) Illustration of a typical adaptive mesh required for the B-F method. The 
fluid-solid interaction applies everywhere in the domain, requiring increased resolution along 
the bottom of the domain, (b) Illustration of the adaptive mesh used for simulations reported 
in this paper. The interface is the only portion of the domain that needs to be resolved at high 
resolution, due to the fact that the fluid-solid interaction is included as an interfacial force. 


C. Long-wave formulation (L-W) 

The long-wave formulation (L-W) describes the evolution of the liquid film thickness, h{x) (or h{x, z) in 
3D), subject to the well-known assumptions of small interfacial slopes (and therefore small contact angles) 
and negligible inertial effects see, e.g. [39] for an extensive review of the L-W. In the present work, we will 
use the linear stability analysis (LSA) of the L-W to validate the simulation results and to highlight the 
differences between the LSA predictions and the direct numerical simulations, in particular in the presence 
of inertia and large contact angles. We note that the focus in this paper is not on extensive comparison 
between the simulation and the L-W results. Such a comparison can be found in [50] for spreading and 
retracting droplets. 

In L-W, an additional contribution, of the form Il{h) = ]CisF{h)^ to the evolution equation is made to 
include the fluid-solid interaction. Note that Il{h) is of the same form as the fluid-solid interaction term, 
Eq. (3), formulated for a flat film. See [41] for more details, and [39, 51] for reviews; an extensive discussion 
regarding inclusion of disjoining pressure in the L-W equation can be found in [18]. With Il{h) included, the 
dimensionless L-W equation is given by 

3ht + Vx,z ■ {h^V ■ [h^Vx,zKF{h)\ = 0 , ( 11 ) 

where Vx,z stands for 2D in-plane gradient operator. Equation (11) is derived using the same scales as the 
ones used to obtain Eq. (8), with the same length scale chosen for both the in-plane directions {x and z) 
and the film thickness {h). The LSA of Eq. (11) shows that if the initial film thickness, is perturbed by 
a mode with infinitesimal amplitude, (5, of the form exp(i(/cx + / 2 ))), then the initial perturbation will grow 
or decay with a growth rate 


hi{e + p){ki-{e + p)) 


( 12 ) 


where kc is the critical wavenumber, given by 


K 


ho 


m 


h*^rn 

— n 


ho 


ho 


(13) 


Note that if < 0, then /3 < 0, and there is no instability for any wavenumber. If > 0, all modes with 
wavenumber k < kc are unstable. Associated with Eq. (12) is a wavenumber of maximum growth, fcmax) and 
the corresponding maximum growth rate, /?niax, given by 


kn 


= ^; /3max=/3(fc' + ^" = fcLx) = §- 


(14) 


We will frequently make reference to the wavelength of maximum growth, defined by Amax = 27r//cn 
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A 

R-P 

B-F 

1/2® 

6.1 X 10-2 

7.0 X 10-3 

1/2® 

2.6 X 10-3 

4.3 X 10-3 

1/2^ 

8.0 X 10-^ 

3.8 X 10-^ 


TABLE I: Comparison of the mesh convergence between the R-P method and the B-F 
method of [41] for a drop. Here, A is the grid resolution used (the maximum one for the R-P 
method, and uniform for the B-F method), and the norm of the error in the profile shape. 
The performance of the two methods is comparable at higher resolutions. 


h* 

R-P 

B-F 

0.03 

1.36 

1.37 

0.015 

1.45 

1.47 

0.0075 

1.49 

1.53 


TABLE II: The contact angle, ^num, obtained by the R-P and the B-F methods for different 
values of h*. The contact angles approach ^eq = 7r/2 at a similar rate. 


Within the L-W, these results imply the following features of the film instability: 

• The wavenumber of maximum growth, /Cmax, scales with yT^^cos^^. Large contact angles imply 
larger /Cmax, and a corresponding decrease in Amax- 

• The maximum growth rate, /^max, scales with (1 — cos 6>eq)^. Large contact angles dramatically increase 
the growth rate of the dominant mode, and thus reduce the time it takes for films to break up. 


III. RESULTS 
A. Contact Angles 

We first demonstrate that the method under consideration can yield contact angles as accurately as the 
methods presented in [41]. For this purpose, we consider a 2D droplet, with the initial fraction (the initial 
region occupied by the liquid phase, i.e. where T = I) set to be the following 

{(x, y) \ ^ {y ^ RcosOi — < B? or y < /i*}. 

Thus the initial profile is a circular cap sitting on top of an equilibrium film such that it intersects this film 
with angle Oi. The value of R is chosen so that the total area of the circular cap is equal to Aq = 0.75‘^7rl2. 
For simplicity, for our first tests, we set = Oi = 7r/2. Since for small /i*, 6>num is close to 6>eq, the initial 
fluid configuration is close to its equilibrium. The droplet is then allowed to relax, and we measure ^num via 
the method described in Sec. IT We set Oh = (0.05)“^/^, a sufficiently large value so that inertial effects are 
not significant. 

Table I compares the convergence of the equilibrium droplet profile as a function of the minimum mesh 
size, for a droplet simulated using the R-P method discussed here, and the B-F method of [41]. The R-P 
simulations are resolved to a resolution Amax = A on an adaptive mesh as shown in Fig. I, while the B-F 
method simulations are refined uniformly with mesh size A. The error is measured as the L^ norm of the 
error in the profile shape. The equilibrium film thickness is /i* = 0.03. Despite the fact that the B-F 
simulations are run with a uniform mesh, the R-P method performs comparably well in convergence in mesh 
to the B-F method at higher resolutions. 

Figure 2 shows simulation profiles obtained using the R-P method, again with 6 >eq = Oi = 7 r/ 2 . The 
profiles at equilibrium are shown for /i* = 0.03 (red) , /i* = 0.015 (green), and /i* = 0.0075 (blue). The 
initial condition for /i* = 0.015 is plotted by the black dotted line, showing the profile of a droplet with 
contact angle ^eq- The smooth transition from the droplet profile to the equilibrium film is clearly visible; 
as /i* is reduced, the equilibrium profiles approach those of a droplet with contact angle O^q. 
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FIG. 2: Effect of h* on the profile when Oi = ^eq = 7r/2. The dotted line shows the initial 

profile for h* = 0.015. 


Table II shows the measured values of 6>num, agaiu for the droplet with 6>eq = Oi = 1^12^ as /i* is varied, 
for both the B-F aud the R-P methods. While the values differ slightly betweeu the two methods, as /i* is 
reduced, 6>num becomes closer to 6>eq. 

We couclude that the R-P method is comparable to the B-F method iu its ability to model the coutact 
augle. Iu additiou, it comes with the poteutial for dramatically improved performauce, due to fact that the 
layer of fluid uear the substrate does uot uecessarily ueed to be highly resolved. For example, for the drops 
simulated iu Table I, with the timestep fixed at 10“^, at A = 1/2^ the R-P method takes approximately 14% 
as much CPU time as the B-F method; at A = 1/2^, the R-P method has approximately 5% the ruutime. 
This is because, for each timestep, the uumber of computatious each method iucurs is 0{N), where N is the 
uumber of cells. However, the adaptive mesh illustrated iu Fig. 1 has approximately 0(1/A) cells, while a 
uuiform mesh has 0(1/A^) cells. Simple adaptive meshes were used iu [41], improviug the performauce of the 
B-F method, however, higher resolutiou is still required at a layer of uouzero thickuess uear the substrate. 
The complexity differeuce is similar iu 3D, where the R-P method scales as 0(1/A^) aud the B-F method 
as 0(1/A^). 

The improvemeut iu performauce permits the study of problems that are impractical to cousider usiug the 
B-F method, iucludiug breakup of the films iu 2D aud 3D. We proceed to aualyze these problems. 


B. Film Instability: Linear regime 

Iu this sectiou, we compare the LSA predictious of the L-W equatiou with the results of simulatious, both 
applied to the iustability of thiu films iu 2D. We focus here ou early times/liuear regime for which the LSA 
is expected to hold aud discuss the comparisou betweeu the results iuflueuced by the relevaut parameters. 
Iu particular we focus ou the iuflueuce of the uuperturbed film thickuess, ho, aud of the equilibrium coutact 
augle, ^eq- Note that the LSA should apply eveu for large coutact augles, siuce we focus ou early stages 
of iustability wheu the iuterfacial slopes are small. We set Oh = 0.45, which implies a low iuertia regime, 
so that the assumptious of the L-W are reasouably well satisfied (iu Sec. IIIC, we discuss the values of Oh 
that coustitute low iuertia iu the preseut coutext). For all the simulatious, the iuitial regiou occupied by the 
liquid phase (i.e., T = 1) is betweeu ^ = 0 aud h = h{t = 0,x) = ho{l + Scos{kx)), where 6 = 0.08. For all 
simulatious iu this sectiou, we set /i* = 0.12. 

Figure 3 shows the growth rates resultiug from the LSA together with the oues extracted from simulatious 
for early time evolutiou, for a discrete set of waveuumbers. We measure the growth rate of simulatious 
accordiug to the followiug procedure: first, we output the miuimum ^-coordiuate of the fluid iuterface, 
hmin{t)' We theu fiud to aud U, such that for to <t<ti, the fuuctiou \og{hmin{t)) is approximately liuear 
iu t; to is geuerally uear 0, aud ti is small, so this correspouds to the iuterval of time wheu the growth of the 
uustable modes is goverued by the LSA. We theu perform a least-squares fit of a Hue to log{hmin{t)) over 
this iuterval. Repeatiug the same procedure for the maximum ^-coordiuate of the fluid iuterface, hmax{t), 
we calculate the growth rate as the average of the slopes of the fitted Hues. 

We focus ou tHe iuflueuce of ho, aud set 6>eq = 7r/2. For small waveuumbers, k, the agreemeut is very good; 
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FIG. 3: Comparison between the dispersion curve predicted by the LSA of the L-W equation 
(blue solid curve) and the growth rate computed by R-P based simulations (symbols), for 
^eq = 7r/2, for (a) ho = 1.0 and (b) ho = 0.5. 


however, for larger values of h, there are significant differences, for both values of ho- We conjecture that 
the differences are due to the fact that the L-W assumptions no longer hold for larger values of h, since the 
separation of scales required for the L-W to apply may not be satisfied. For example, if one considers the 
relevant lengthscales as Xc = Amax and he = ho ^ in the in-plane and out-of-the-plane direction, respectively, 
then for ho = 0.5 one finds e = Xc/hc is 0(0.1) and therefore for k values larger than hmax, the requirement 
e ^ 1 is no longer a good approximation. However, it should be noted that the differences between the 
simulations and the LSA are mostly focused on large values of h, where films are stable. 

One may expect that the LSA would lead to more accurate results for thinner films. As we can see in 
Fig. 3, this is not the case. An insight can be reached again by considering the relevant length scales; the 
ratio of the film thickness, ho, and the relevant in-plane length scale, Amax, given by 


he 

Amax 



(15) 


This ratio is a decreasing function of ho for ho > h* [m{m — l)/{n{n — 1))]’^”’^. Therefore, reducing ho 
actually implies that the LSA is less accurate, except when ho is nearly as small as h*. 

Equation (15) implies that a reduction in ^eq should improve agreement with the predictions of the LSA. 
Figure 4 shows the results for three values of ^eq and we indeed immediately observe much smaller differences 
between the simulations and the LSA for smaller 6>eq. For larger 6>eq, the LSA strongly overpredicts the growth 
rates as well as the values of hmax- We will see in the next section that the differences between the values 
of ^max resulting from the simulations and the LSA become much more visible in the nonlinear regime of 
instability development. 


C. Film instability: Nonlinear evolution and breakup - 2D 

Here we focus on the nonlinear process of film breakup, and concentrate in particular on the role of the 
contact angle, 6>eq, and Ohnesorge number, Oh, on the properties of emerging patterns. While the discussion 
that follows is general, and is formulated in terms of previously defined nondimensional quantities, in order to 
choose a reference parameter set, we consider nanoscale liquid metallic films. The stability properties of such 
films have been a subject of recent interest in connection with nanoparticle fabrication (see e.g. [1, 3, 30-33]). 
In experiments, a film of metal, with thickness on the order of tens of nanometers, is deposited on a surface. 
This film is then liquefied, and subsequently breaks up into droplets, which solidify to form nanoparticles. 

For definitiveness, we focus on Cu films, as considered in [4]. The viscosity, density, and surface tension 
are taken to be the values of liquid Cu at its melting point [4], with p = 7760 kg/m^, p = 0.00438 Pa-s, 
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FIG. 4: Comparison between the growth rate predicted by the LSA (blue solid curve), and the 
simulations (symbols), for ho = 0.5, for (a) ^eq = 37r/4, (b) ^eq = 7r/2, and (c) ^eq = tt/G. The 
value of hmax obtained from simulations is shown by the vertical solid line; the dash-dotted 
line shows the value of hmax predicted by the LSA. 


and 7 = 1.304 N/m. We set relevant length scale, L, as the reference film thickness, taken to be 8 nm. 
The Ohnesorge number corresponding to these parameters is Oh = 0.487. Note that this represents a rough 
lower bound for Oh in the motivating experiments: temperatures may exceed the melting point substantially, 
reducing the viscosity, and furthermore the relevant length scales may be larger than 8 nm. This variation in 
Oh makes it important to explore the boundary between the viscous and inertial regimes. The simulations 
that we present in this section are for a film of thickness ho = 1; we also considered films with ho = 0.5, 
and obtained similar results. We use /i* = 0.1225, and 6 >eq = 79° ^ 0 . 4397 r. Note that for these parameters, 
Amax = 15.6, and Ac = 11.0. 

For 2D simulations, the computational domain is specified by (Lx^Ly), with Lx ~ 8 Aniax 5 and Ly suffi¬ 
ciently large so that the value assigned to it does not influence the results. For the simulations discussed 
next, we use Lx = 122.5, Ly = 43.125, and the initial film shape specified by is ho + C (^)5 where 

C(a;) = ^<5iCos (^y (16) 

i=l \ W 

Here, A^ = 2Lxli^ and hi is a random perturbation amplitude, uniformly distributed in the range ±0.0125. 
We simulate the evolution with Ns = 20 sets of hi. For each simulation (numbered j), a discrete height 
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profile is produced, hj{t,x). We then compute the discrete Fourier transform (DFT) of each height profile, 
Hj{t, k). Finally, we compute the average of these as 


-j s 

H2D{t,k) = ^Y.^iXk)- (17) 

Figure 5 shows the time evolution of a typical profile hi{t^x) (the top image in each of Fig. 5(a) - (f)), 
H 2 D{t^k) along with the dispersion curve from the LSA (the bottom image in each of Fig. 5(a) - (f), 
smoothed with a running 5 point average). We see that as the initial perturbations begin to grow, the 
H 2 D{t,k) attains the form similar to the dispersion curve (Fig. 5(a) - (c)). However, as the holes and 

consequently drops start to form, the peak in H 2 D{t^ k) shifts towards smaller values of k, so that the final 
distribution of drops is characterized by a length scale larger than Amax- 
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FIG. 5: Time evolution of a representative simulation of film breakup in 2D with ho = 1, 
Oh = 0.487, and ^eq = 0.4397r; the solid blue curve shows the fluid interface. The associated 
Fourier spectrum, averaged over 20 realizations, is shown below each image, (a) t = 0, (b) 
t = 179, (c) t = 357, (d) t = 536, (e) t — 715, and (f) t = 882. 
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FIG. 6: Time evolution of a representative simulation of film breakup in 2D with ho = 1, 
Oh = 0.0487, and ^eq = 0.4397r; the solid blue curve shows the fluid interface. The associated 
Fourier spectrum, averaged over 20 realizations, is shown below each image, (a) t = 0, (b) 
t = 357, (c) t = 536, (d) t = 7l|3 (e) t = 1429, and (f) t = 2858. 







































































FIG. 7: Comparison of as a function of time for films of varying Oh. Symbols show the 
approximate breakup times for each parameter set, and the solid dashed line shows /cmax 
predicted by the LSA. Note that Oh = 4.87 and Oh = 0.487 are visually nearly 
indistinguishable. The ticks on the y-axis are for wavenumbers that can be resolved on the 
finite domain. Here, /lo = 1, ^eq = 0.4397r. 


Next we consider different values of Oh. For Oh = 4.87, no significant difference is observed compared 
to Oh = 0.487, indicating that both sets of results belong to the high viscosity limit. Smaller values of Oh, 
however, lead to different results. Figure 6 shows profiles and associated averaged Fourier spectra for a film 
with Oh = 0.0487. We see that the film takes much longer to break up compared to the films characterized 
by larger Oh. Also, the long term evolution of the Fourier spectrum shows a flatter peak, when compared 
with the larger Oh simulations, indicating that the preference for a specific wavenumber is weaker in the late 
stages of evolution. The practical consequence of this may be increased degree of disorder in the distribution 
of emerging patterns. 

In order to study the effect of Oh on the evolution of the film more systematically, we define as 

the value of k for which H 2 D{k) attains its maximum. Figure 7 plots k^^ for various Oh as a function 
of time (note that due to the finite domain size, k^^ only takes discrete values, leading to the staircase 
appearance of these plots). The symbols on each curve indicate the approximate time at which the film 
ruptures, and the dashed line shows the value of kma^ from the LSA. The overall behavior is similar for all 
curves: k^^ quickly takes a value as close to the /cmax as can be resolved on a finite domain, and then stays 
at this vaiue untii the him breaks up. Afterwards, it decreases, indicating that the iength scaie of the finai 
dropiet distribution is iarger than Amax- We note that for Oh = 4.87 and Oh = 0.487, the behavior is neariy 
indistinguishabie, indicating large Oh limit. Most importantly, smaller values of Oh significantly increase 
the time it takes for the film to rupture. The LSA does not capture this behavior, and instead predicts that 
the breakup time is approximately the same (~ 500) for all values of Oh, consistently with the results of 
large Oh simulations. 

Next, we investigate the effect of the contact angle on the rupture. We carry out two additional simulation 
sets, both with Oh = 0.487, and with different values of 6>eq and different domain sizes, keeping Lx ~ SAmax- 

1. 6>eq = 7r/6, Lx = 367.5; 

2. 6>eq = 37r/4, Lx = 91.875. 

For each set, the form of the perturbation is the same as in Eq. (16), using corresponding values of Lx. 
Figure 8 shows k^^ for these two contact angles, as well as for ^eq = 0.4397r. For each curve, the symbol 
of the same color indicates the breakup time, and the dashed line of the same color shows /Cmax from the 
LSA. Note that the breakup times predicted by the LSA are « 18000, 500, and 110, for ^eq = 7r/6, 0.4397r, 


14 


























FIG. 8: Plot of for varying Oeq- The green curve corresponds to Oh = 0.487 in Fig. 7. 
The symbols show the approximate time of him rupture. Each dashed line shows /cmax from 
the LSA for the curve of the same color. 


and 37r/4, respectively, roughly in line with the simulated breakup time for all cases. For 6>eq = 37r/4, the 
behavior is similar to ^eq = 0.4397r. However, the evolution of ^eq = 7r/6 is dramatically simpler: almost 
immediately, relaxes to near /cmax, and remains at approximately the same value for the considered 

simulation times. 

To gain a basic understanding of the influence of ^eq on the evolution, consider the following simple 
argument. Approximate the prohle of a rupturing him of unit initial thickness, perturbed by a mode of 
wavelength Amax, by h{x) = 1 + A{t) cos (27rx/Amax), where A{t) is the time dependent amplitude. At the 
time of breakup, A ^ 1. Then, further approximate the instantaneous contact angle. Or, of the rupturing 
him by considering the slope at the point of inhection, given by tan^^ = 27r/Amax- Substituting the value of 

Amax, we obtain tan^^ = By/1 — cos^eq, with B = 2\/2/i*. Consider now the difference 

between Or and 6>eq. This difference is an increasing function of 6>eq, suggesting that when hlms with larger 
6>eq rupture, the corresponding Or is signihcantly different from 6>eq, leading to quick retraction of the him 
from the point of rupture. This retraction reduces the space available for consecutive breakups, and leads 
to larger distances between the drops than predicted by the LSA. For smaller values of ^eq, the difference 
between Or and 6>eq is small at the point of initial rupture, and the retraction effect is not present. 

Next we comment on the dihFerence between k'^^ and the corresponding kmaxi visible in Figs. 7-8. First 
we note that this shift is not affected by inertial effects for the simulated times: Fig. 7 shows that Oh = 48.7 
and Oh = 0.0487 lead to the same k'^^ after breakup (the additional coarsening expected for very long 
times is beyond the scope of the present study). These results cannot be explained by the breakdown of 
the LSA for large Oeq’ from Fig. 4, we might expect a roughly 10% difference between the wavelength of 
maximum growth in LSA and simulations, which is significantly less than the difference between k'^^ and 
kmax after the film has ruptured. Thus, we expect that the difference is due to nonlinear effects relevant 
during breakup for larger contact angles. We note here that for for smaller contact angles the effect is not 
visible within the present degree of accuracy; the simulations carried out in larger domains and with larger 
number of realizations using long-wave approach suggest that (much less obvious) coarsening effects may be 
seen within long-wave simulations as well [52]. 
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FIG. 9: Time evolution of a representative simulation of film breakup in 3D, with Oh = 0.487, 
ho = 1.0, and Oeq = 0.4397r. The color shows the logarithm of the height of the interface above 
the substrate. The associated Fourier spectrum is shown below each image; these data are 
averaged over 10 instances, and smoothed with a 5 point running average, (a) t = 0.0, (b) 
t = 335, (c) t = 558, and (d) t = 781. 
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FIG. 10: Comparison of as a function of time for 3D films (blue) and 2D films (green). 
The green inverted triangle shows the approximate breakup time of the 2D simulation. The 
blue triangle shows the approximate time at which the first holes form in the 3D film; the 
inverted blue triangle shows the approximate time that droplets begin to form. The solid 
dashed line shows /cmax predicted by the LSA. 


D. Film instability: Nonlinear evolution and breakup - 3D 


Next, we consider film breakup in 3D. The setup is identical to the 2D case, except that the domain size 
is specified by = 4Aniax and Ly = 31.25. The initial condition is a perturbed film specified by 

ho + C(x, z), and 


i=l j=l ^ ^ ^ \ / 

where and are random perturbation amplitudes, uniformly distributed in the range ±0.0125. 

We simulate the film for Ns = 10 sets of As in the 2D case, we produce a height profile for each 
simulation, hj{t^x^z)^ compute its discrete Fourier transform, and the average 


Hsnit, k) 


N. 


Ns 




i=i 


(19) 


The summand in Eq. (19) represents the average of the discrete Fourier transforms along the coordinate 
axes, taking advantage of the symmetry of Eq. (12) to get more information about the growth rates. 

Eigure 9 shows the time evolution of hi{t^x^z) (top image in each part of the figure), and H^d along 
with the dispersion curve predicted by Eq. (12) (bottom images). The plots of Hod have been smoothed 
by a 5 point running average. A similar trend is observed as in Eig. 5: the spectrum takes on a similar 
profile to the dispersion curve from the LSA, and as breakup proceeds, its peak, k^^, shifts towards smaller 
wavenumbers. 

Eigure 10 shows k^^{t) for both 2D and 3D films (Oh = 4.87). Eor both 2D and 3D, k^^ is close to 
hmax until the film ruptures, after which both k^^ relax to smaller values. However, once droplets begin to 
form in 3D simulations, k^^ relaxes to an even smaller value. The same retraction argument that we used 
to explain the evolution of k^^^(t) in 2D can as well be used here to explain this phenomenon. 

We also study the effect of initial film thickness on the evolution of the Eourier spectra for 3D films. 
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Figure 11 presents the evolution of the film with an initial thickness ho = 0.5, while keeping Oeq, h*, and Oh 
the same as in Fig. 9. Our initial condition is again subject to a white noise perturbation (Fig. 11(a)). Holes 
form rapidly (Fig. 11(b)), due to the larger growth rate associated with the smaller ho, when compared with 
Fig. 9. Hsd shows a peak near kmax at the time of the formation of holes; however, this does not agree as 
closely with the overall shape of the dispersion curve as is observed for thicker films. We conjecture that this 
difference is due to the fact that holes form almost immediately, so that there is insufficient time for Hsd to 
relax to the shape of the dispersion curve. Similarly to thicker films, as drops begin to form, relaxes 

to wavenumbers significantly smaller than kmax- 

Fig. 12 shows the time evolution of k'^^ for 3D films with Hq = 1 and ho = 0.5. Similarly to the case 
when ho = I plotted in Fig. 10, k^^^ shifts to smaller wavenumbers after two distinct events: first, when 
holes begin to form in the film, and second when drops begin to form (marked by triangles and inverted 
triangles, respectively). However, two significant differences are observed for the thinner film. First, breakup 
takes place almost immediately, so there is no interval of time where k'^^ approximates kmax- Second, for 
larger times, the difference between kmax and k^^^ is greater for ho =0.5 than for ho = 1. 

To summarize, the DFT of the nonlinear film breakup shows that the dominant length scales deviate from 
Amax predicted by the LSA. The degree of deviation depends on heq, and whether the film is 2D or 3D; 
additionally, the time it takes for the film to rupture depends strongly on Oh. For 2D simulations with 
small heq, the DFT has a peak at or close to hmax for the entirety of the film evolution. For larger heq, the 
evolution of the DFT shows two distinct phases: prior to breakup, when the peak in the DFT is at k^ax^ and 
after breakup, when the peak shifts to smaller wavenumbers. For 3D simulations, the evolution of the DFT 
shows three distinct phases, each associated with a shift towards smaller wavenumbers: prior to breakup, 
after holes begin to form, and after drops begin to form. A decrease in Oh is primarily associated with a 
dramatic increase in the time it takes for the film to rupture. 
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FIG. 11: Time evolution of a film breakup in 3D, with /lo = 0.5 and all parameters otherwise 
identical to the simulation presented in Fig. 9. The color shows the logarithm of the height of 
the interface above the substrate. The associated Fourier spectrum is shown below each 
image; these data are averaged over 10 instances, and smoothed with a 5 point running 
average: (a) t = 0.0, (b) t = 45, (c) t = 134, and (d) t = 301. 


IV. CONCLUSIONS 


In this paper, we have demonstrated a computationally efficient method for including fluid-solid interac¬ 
tions into direct numerical simulations. This method is found to perform as well as the body force formula¬ 
tion [41], while requiring a significantly smaller computational effort. The two methods were compared by 
considering contact angles of equilibrium drops in 2D, where both methods perform similarly in terms of 
convergence with mesh refinement, as well as in terms of the behavior when the equilibrium film thickness, 
/i*, is reduced. The computational complexity required is however dramatically reduced for the presently 
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FIG. 12: Comparison of as a function of time for 3D films with initial thickness Kq — 1 
(blue) and ho = 0.5 (green). The blue curve is identical to that seen in Fig. 10. The blue 
(green) triangles show the times at which the first holes appear, and the inverted ones show 
the times at which droplets begin to form. The blue (green) dashed lines show the 
corresponding values of hmax- Oh = 0.487, and ^eq = 0.4397r for both simulation sets. 


considered method. 

With the established improvement in computational performance, we are now able to study the instability 
of films due to fluid-solid interaction using direct numerical simulations to gain a quantitative understanding 
of the evolution of the instability, film rupture, and the post-rupture dewetting process. We compare the 
results of direct simulations with the linear stability analysis (LSA) of the long-wave formulation (L-W), and 
find that when contact angle, is larger, there is a significant difference with the predictions of the LSA. 
We also demonstrate that a reduction in the film thickness does not reduce the differences between the LSA 
and direct simulations, as the typical in-plane lengthscale decreases even faster than the out-of-plane one 
(film thickness). 

Finally, we study breakup of a film in both 2D and 3D. We consider the evolution of the emerging length 
scales by computing the discrete Fourier transform (DFT) of the fluid profile. 2D films are characterized 
by a two stage evolution, and 3D films by a three stage one; for both cases, the initial phase exhibits a 
DFT with a peak, near /cmax from the LSA, and each successive stage is associated with a decrease 

in and correspondingly, an increase in the emerging length scales. When a flat film is perturbed by 

small perturbations, the simulations illustrate the dynamics of dewetting, the formation of ‘dry spots’, and 
the subsequent coarsening due to inertial effects. When studying the effect of the contact angle for 2D films, 
we find that the dewetting morphology depends on the contact angle; particularly, the agreement between 
the simulation results and the predicted fastest growing unstable mode from the LSA deteriorates for large 
contact angles after the formation of the first expanding holes. When studying the breakup of films in the 
presence of significant inertial effects, so that characteristic Ohnesorge number. Oh, is small, we find that 
the long term preference for a specific wavenumber is weaker when compared with the larger Oh, resulting 
in an increased degree of disorder in the distribution of emerging patterns. 

The speed and simplicity of the implementation of the method presented in this paper opens up the 
possibility of studying a variety of problems involving thin films and other geometries that have not been 
studied extensively so far by direct numerical simulations. Our numerical method allows, for the first time, 
for the simulation of arbitrary contact angles (including those greater than 7r/2), and film breakup, including 
inertial effects. While we have studied the influence of varying contact angle and Oh in 2D films, we leave 
the exhaustive parameter study of 3D films for a future work. 
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